/* Copyright (c) Colorado School of Mines, 2010.*/
/* All rights reserved.                       */

/* SUEA2DF: $Revision: 1.26 $ ; $Date: 2010/02/04 18:44:43 $		 */

#include "su.h"
#include "segy.h"
#include "header.h"

/*********************** self documentation **********************/
char *sdoc[] = {
"									",
" SUEA2DF - SU version of (an)elastic anisotropic 2D finite difference 	",
"		forward modeling, 4th order in space			",
"									",
" suea2df > outfile c11file= c55file  [optional parameters]		",
"									",
" Required Parameters:							",
" c11file=c11_file	c11 voigt elasticity parameter filename		",
" c55file=c55_file	c55 voigt elasticity parameter filename		",
"									",
" Optional Parameters:							",
" rhofile=rho_file	density filename				",
"			(if rhofile is not set, rho=1000 is assumed)	",
" Anisotropy parameters:						",
" aniso=0	 	 =1 - include anisotropy parameters		",
" mode=0		=0 output particle velocity, =1 output stresses ",
"			(snapshots only)				",
"									",
" ... the next 3 parameters become active only when aniso=1....		",
" c13file=c13_file	c13 voigt elasticity parameter filename		",
" c33file=c33_file	c33 voigt elasticity parameter filename		",
" c15file=c15_file	c15 voigt elasticity parameter filename		",
" c35file=c35_file	c35 voigt elasticity parameter filename		",
"									",
" Attenuation parameters:						 ",
" qsw=0		 switch to include attenuation =1 - include		",
" ... the next parameter becomes active only when qsw=1....	     	",
" qfile=Q_file	  Q parameter filename	    				",
"									",
" dt=0.001		time sampling interval (s)			",
" ft=0.0 		first time (s)				 	",
" lt=1.0 		last time (s)					",
"									",
" nx=200		number of values in slow (x-direction)		",
" dx=10.0	 	spatial sampling interval (m) x-coor		",
" fx=-1000		first x coor (m)				",
"									",
" nz=100		number of values in fast (z)-dimension		",
" dz=dx			spatial sampling interval (m) z-coor		",
" fz=0			firstz coor (m)  				",
"									",
" Source parameters:							",
" sx=0			source x position (m)				",
" sz=500		source location (m)  				",
" stype='p'		source type					",
"			  p: P-wave					",
"			  v: velocity					",
"			 pw: P plane-wave				",
" sang=0		for stype='pw': plane wave angle		",
" wtype='dg'		wave type					",
" 			 dg: Gaussian derivative 			",
" 			 ga: Gaussian		 			",
" 			 ri: Ricker					",
" 			 sp: spike, sp2: double spike   		",
" ts=0.05		source duration (s)				",
" favg=50		source average frequency			",
"									",
" Attenuation parameters:						",
" qsw=0		 	switch to include attenuation =1 - include	",
"									",
" Boundary condition parameters:					",
" bc=10,10,10,10 	Top,left,bottom,right boundary condition	",
" 			=0 none						",
" 			=1 symmetry 					",
" 			=2 free surface (top only)			",
" 			>2 absorbing (value indicates width of absorbing",
"			layer	 					",
" bc_a=0.95;		bc initial taper value for absorbing boundary  ",
" bc_r=0.;		bc exponential factor for absorbing boundary  	",
" 			variables are scaled by bc_a*pow(i,-bc_r)	",
"									",
" Optional output parameters:						",
" sofile=		name of source file				",
" snfile=		name of file containing for snapshots		",
" snaptime=		times of snapshots i.e. snaptime=0.1,0.2,0.3	",
"									",
" vsx=			x coordinate of vertical line of seismograms	",
" hsz=			z coordinate of horizontal line of seismograms	",
" vrslfile=\"vsp.su\"	output file for vertical line of seismograms[nz][nt]",
" hsfile=\"hs.su\" 	output file for horizontal line of seismograms[nx][nt]",
" tsw=0		 	switch to use shear stress only in non-fluid	",
"			media - may help reduce dispersion tsw=1. If	",
"			tsw=0 then standard calculation	  		",
" verbose=0		=1 to print progress on screen			",
"									",
" Notes:								",
" 1) The outfile contains information generated by the input parameters,",
"    such as memory allocation, stability, etc. If your input file does	",
"    not work, check this file first.					",
"									",
" The model is specified as binary files of stiffness parameters and    ",
" densities. These may be created any way the user desires. The program ",
" unif2 or makevel may be used to generate densities, and the program	",
" unif2aniso may be used to generate the stiffnesses. You will need to	",
" transpose these files (stiffnesses and densities), as the input	",
" format for suea2df assumes that the fast dimesion is the horizontal or ",
" the x-dimension. You may do this via					",
" 									",
"  transp n1=nz < c11_file > transp_c11_file				",
"									",
" If aniso=1 then the program will expect the additional stiffnesss files.",
"									",
" If qsw=1 unif2anis can be used to generate the Q values on a grid	",
" These value also need to be transposed, as with the stiffnesses.	",
"									",
"Output files (always generated)					",
"	hsfile								",
"	vrslfile							",
"	hsfile.chd - header for hsfile					",
"	vrslfile.chd - header for vrslfile				",
"	hsfile.mod - model file						",
"									",
" Output files (if requested)						",
"	sofile - ascii source file					",
"	snfile  - su format snapshots file				",
"									",
" Caveat:								",
" A common error in using this program is to compute stiffnesses with	",
" a specified density, but forget to specify this density as the rhofile.",
NULL};

/* 
 * Credits: UU GEOPHY Chris Juhlin 15 May 1999
 * Copyright (c) Uppsala University, 1998.
 * All rights reserved.			
 * Parts of program use Seismic Unix Package - CSM
 * Changes - C. Juhlin
 *
 * 1. Fixed upgrading of stresses. There was an error in the coding for
 * the Tzz term, c15 was being used instead of c35. This only caused
 * problems for dipping anisotropic layers
 *
 * 2. Added some header information for hutput of snapshots.
 *
 * 3. 2001-01-30: Added option to set absorbing bc constants bc_a and bc_r 
 *
 * 4. 2001-02-23: Corrected bug in outputting model boundaries to standard
 * output in 
 * routine get_econst
 *
 * 5. 2001-04-26: Added option for updating velocities to only use 
 * shear stress if material is non-fluid, this appears to reduce dispersion at 
 * near grazing angles for fluid-solid boundary. Set tsw=1 to invoke
 *
 * 6. 2001-05-14: Changed loop in free-surface boundary condition for velocties
 * Thanks goes to Mike Holzrichter for pointing out this problem and the wrong
 * scaling factor in the updating.
 *
 * 7. 2001-05-16: Changed set_layers function to avoid negative indexing.
 * Thanks goes to Mike Holzrichter for pointing out this problem
 *
 * 8. 2001-05-17: Modified make_seis to take into account VSP geometry
 * correctly and not store unnecessary data.
 *
 * 9. 2001-08-21: Fixed set_layers so mode fills properly in depth. Earlier
 * versions were accessing incorrect array locations at last defined depth.
 *
 * 10. 2003-04-21: Fixed boundary conditions.
 *
 * 11. 2003-05-02: Extended the model area by half the grid spacing on the RHS.
 * This makes the model area symmetric allowing a plane wave source to be
 * introduced into the model (stype=pw). The w, txx and tzz grids contain now
 * one more column than the u and txz grids.
 *
 * 12. 2003-05-02: Added routines to allow plane waves to be introduced at a 
 * specified angle (sang) into the model with functions add_pw_source_V and
 * add_pw_source_S.
 *
 * 13. 15 Oct 2005 -- tossed all the model building stuff. Read models
 *	from binary files made by  unif2aniso (CWP:John Stockwell)
 *
 * 14. 25 Feb 2008 -- Fixed attenutation option (qsw=1) so that Q values are
 *	from binary files made by makevel or similar program
 *
 * Algorithm based on Juhlin (1995, Geophys. Prosp.)
 *	and Levander (1988, Geophysics)
 * Attenuation included as in Graves (1996, BSSA)
 *
 */

/**************** end self doc ********************************/

#define	F1	1.125
#define	F2	-0.04166667

segy tr, trv, trh;

/* prototypes for functions defined and used below */
int get_source(float dt, float ts, float favg, char *wtype, float *source);
void update_vel(int ifx, int ilx, int jfz, int jlz,
	float **u, float **w, float **txx, float **tzz, float **txz,
	float **rho, float **c55, float dtx, float dtz);
void update_vel_tsw(int ifx, int ilx, int jfz, int jlz,
	float **u, float **w, float **txx, float **tzz, float **txz,
	float **rho, float **c55, float dtx, float dtz);
void update_stress_iso(int ifx, int ilx, int jfz, int jlz,
	float **u, float **w, float **txx, float **tzz, float **txz,
	float **c11, float **c55, float dtx, float dtz);
void update_stress_ani(int ifx, int ilx, int jfz, int jlz,
	float **u, float **w, float **txx, float **tzz, float **txz,
	float **c11, float **c55, float **c33, float **c13, 
	float **c15, float **c35, float dtx, float dtz);
void add_p_source(float **txx, float** tzz, float amp, int i, int j);

void add_v_source(float **u, float **w, float amp, int i, int j);

void add_pw_source_V(float **u, float **w, float **txx, 
			float **tzz, float **txz, float **c11, 
			float **c55, float *a, float sang, 
			int il, int ir, int jt, int jb, 
			int ns, int k, float dx, float dz,
			float dt, float *vl, float *vb, 
			float *vr, float *ttl, float *ttb, float *ttr);
void add_pw_source_S(float **u, float **w, float **txx, 
			float **tzz, float **txz, float **c11,
			float **c55, float *a, float sang,
			int il, int ir, int jt, int jb,
			int ns, int k, float dx, float dz,
			float dt, float *vl, float *vb,
			float *vr, float *ttl, float *ttb, float *ttr);
void make_snap(float **u, float **w, float sx, float sz, float dx, float dz,
			int nxpadded, int nzpadded, int nx, int nz, int k,
			float t, float fx, float fz, segy tr, FILE *sneisfp,
			int *wbc);
void make_seis(float **ut, float **wt, float sx, float sz, float dn, float dt,
			int nn, int nt, float dd, float fx, float fz, segy tr, 
			FILE *fp, int b1, int b2, int ivh, 
			int verbose, float favg);
void abs_bc_top(float **u, float **w, float **txx, float **tzz,
			float **txz, int nxpadded, int nzpadded, 
			float *r, int *bc, int *wbc, int ifx, int ilx);
void abs_bc_left(float **u, float **w, float **txx, float **tzz, 
			float **txz, int nxpadded, int nzpadded, 
			float *r, int *bc, int *wbc, int jfz, int jlz);
void abs_bc_bot(float **u, float **w, float **txx, float **tzz, float **txz,
	  		int nxpadded, int nzpadded, float *r, int *bc,
			int *wbc, int ifx, int ilx);
void abs_bc_right(float **u, float **w, float **txx, float **tzz, float **txz,
			int nxpadded, int nzpadded, float *r, int *bc,
			int *wbc, int jfz, int jlz);
void fs4s_bc_top(float **txx, float **tzz, float **txz, int nxpadded,
			int nzpadded);
void fs4v_bc_top(float **u, float **w, float **c11, float **c55, 
			int nxpadded, int nzpadded);
void sym4s_bc_top(float **txx, float **tzz, float **txz, int nxpadded,
			int nzpadded);
void sym4v_bc_top(float **u, float **w, int nxpadded, int nzpadded);
void sym4s_bc_left(float **txx, float **tzz, float **txz, int nxpadded,
			int nzpadded);
void sym4v_bc_left(float **u, float **w, int nxpadded, int nzpadded);
void sym4s_bc_bot(float **txx, float **tzz, float **txz, int nxpadded,
			int nzpadded);
void sym4v_bc_bot(float **u, float **w, int nxpadded, int nzpadded);
void sym4s_bc_right(float **txx, float **tzz, float **txz, int nxpadded,
			int nzpadded);
void sym4v_bc_right(float **u, float **w, int nxpadded, int nzpadded);

void calc_area(float vmax, float dtx, float dtz, int *limits, int i,
			int j, int k, int nxpadded, int nzpadded);
void write_chd(int nxpadded, int nzpadded, int nt, float dx, float dz, float dt,
			float fx, float xmax, float fz, float zmax,
			float sx, float sz, float favg, float ts, 
			char *wtype, char *stype, int *bc, int qsw, 
			int aniso, float xz, FILE *fp, int comp, char *mfile);
int fwrite_chd(FILE *fp, char *s, int k);
void write_grid(float **u, float **w, float **txx, float **tzz, 
			float **txz, int is, int js);
void make_stress_snap(float **txx, float **tzz, float **txz, float sx,
			float sz, float dx, float dz, int nxpadded,
			int nzpadded, int k, float t, float fx, float fz,
			segy tr, FILE *sneisfp, int *wbc);
void pad_float_2D(int n1, int n2, int *padval, float **in , float **out); 
		
/* the main program */
int
main (int argc, char **argv)
{
	float dt;		/* time sampling interval 	*/
	float ft;		/* first time sample		*/
	float lt;		/* last time sample		*/

	float dx;		/* spatial sampling interval x 	*/
	float dz;		/* spatial sampling interval z	*/
	float fx;		/* min x coor		 	*/

	float xmax;		/* max x coor		 	*/
	float fz;		/* min z coor		 	*/
	float zmax;		/* max z coor		 	*/

	char *stype=NULL;	/* source type			*/
	char *wtype=NULL;	/* wave type			*/
	float sang;		/* plane wave angle (deg)	*/

	float ts;		/* source length (s)		*/
	float favg;		/* source average frequency (Hz)*/
	float sx,sz;		/* current source location (m)	*/

	int nsnap;		/* number of snapshots		*/
	float *snaptime=NULL;	/* snapshot times		*/
	int *isnap=NULL;	/* snapshot samples		*/

	int nt=0;		/* number of time samples 	*/
	int nz=0;		/* number of z samples of model */
	int nx=0;		/* number of x samples of model	*/
	int nzpadded=0;		/* nz padded by bc		*/
	int nxpadded=0;		/* nx padded by bc		*/

	int nbytes;		/* number of bytes allocated	*/
	int ngrids;		/* number of bytes allocated	*/
	int i,j,k,ix,iz;	/* loop counters 		*/

	float **c11=NULL;	/* elastic constant c11(x,z)	*/
	float **c11temp=NULL;	/* ...related temporary array	*/
	float **c55=NULL;	/* elastic constant c55(x,z)	*/
	float **c55temp=NULL;	/* ...related temporary array	*/
	float **rho=NULL;	/* density rho(x,z), initially	*/
				/* ... later rho(x,z)=1/density	*/
	float **rhotemp=NULL;	/* ...related temporary array	*/
	float **q=NULL;		/* Q-factor q(x,z)		*/
	float **qtemp=NULL;     /* ...associated temp file	*/

	float rhomin=0.0;	/* minimum density		*/
	float rhomax=0.0;	/* maximum density		*/
	float c11min=0.0;	/* minimum c11			*/
	float c11max=0.0;	/* maximum c11			*/
	float c55min=0.0;	/* minimum c55			*/
	float c55max=0.0;	/* maximum c55			*/

	float c13min=0.0;	/* minimum c13			*/
	float c13max=0.0;	/* maximum c13			*/
	float c33min=0.0;	/* minimum c33			*/
	float c33max=0.0;	/* maximum c33			*/
	float c15min=0.0;	/* minimum c15			*/
	float c15max=0.0;	/* maximum c15			*/
	float c35min=0.0;	/* minimum c35			*/
	float c35max=0.0;	/* maximum c35			*/

	float qmin=0.0;	 /* minimum q		    */
	float qmax=0.0;	 /* maximum q		    */

	int qsw;		/* switch to include attenuation*/

	int aniso;		/* switch to include anisotropy */
	float **c13=NULL;	/* elastic constant c13(x,z)	*/
	float **c13temp=NULL;	/* ...associated temp file	*/
	float **c33=NULL;	/* elastic constant c33(x,z)	*/
	float **c33temp=NULL;	/* ...associated temp file	*/
	float **c15=NULL;	/* elastic constant c13(x,z)	*/
	float **c15temp=NULL;	/* ...associated temp file	*/
	float **c35=NULL;	/* elastic constant c35(x,z)	*/
	float **c35temp=NULL;	/* ...associated temp file	*/

	float **u=NULL,**w=NULL;	/* particle velocities		*/
	float **txx=NULL,**tzz=NULL;	/* normal stresses		*/
	float **txz=NULL;		/* shear stress			*/
	float **ut=NULL,**wt=NULL;	/* time sections 		*/

	int tsw;		/* switch for shear stress updating */

	float *source=NULL;	/* source waveform		*/
	float pfac;		/* P-wave source scaling factor	*/
	int ns;			/* number of samples in source	*/

	float energy,efac;	/* energy in model		*/
	float dtx,dtz;		/* fd factors			*/
	int ifx,ilx,jfz,jlz;	/* grid area for operators	*/
	int is,js;		/* source location in grid	*/
	float vlim[2];		/* min and max velocity in grid	*/
	int *limits=NULL;	/* calculation area in grid	*/

	int verbose;		/* is verbose? 			*/

	int vs2;		/* depth in samples of horiz rec line */
	int hs1;		/* horiz sample of vert rec line */
	float hsz;		/* z position of horiz receiver line */
	float vsx;		/* x position of vertical receiver line */

	int bc[4];		/* boundary condition flags 	*/
	int wbc[4];		/* boundary width	 	*/
	int nbc;		/* number of values given 	*/
	float *bc0=NULL;	/* taper at top boundary 	*/
	float *bc1=NULL;	/* taper at left boundary 	*/
	float *bc2=NULL;	/* taper at bottom boundary 	*/
	float *bc3=NULL;	/* taper at right boundary 	*/
	float bc_a;		/* bc initial taper value 	*/
	float bc_r;		/* bc exponential factor 	*/
	float *vl=NULL,*vr=NULL,*vb=NULL;/* velocity at boundaries 	*/
	float *ttl=NULL,*ttr=NULL,*ttb=NULL;/* traveltimes at boundaries */

	float t;		/* current time of propagation	*/


	char *sofile="";	/* source filename		*/
	FILE *soeisfp=NULL;	/* ... its file pointer  */

	char *snfile="";	/* snapshot filename		*/
	FILE *sneisfp=NULL;	/* ... its file pointer  */

	char *mfile=NULL;	/* model filename*/
	FILE *mfp=NULL;		/* ... its file pointer   */
	char *vrslfile="";	/* vert receiver seismogram line file  name */
	FILE *vseisfp=NULL;	/* ... its file pointer  */
	char *vrshfile=NULL;	/* vert receiver seismogram header filename*/
	FILE *vchdfp=NULL;	/* ... its file pointer  */

	char *hsfile="";	/* horiz receiver seismogram line filename*/
	FILE *hseisfp=NULL;	/* ... its file pointer  */
	char *hcfile=NULL;	/* horiz receiver seismogram header filename*/
	FILE *hchdfp=NULL;	/* ... its file pointer   */
	char tfile[256];	/* temporary filename*/

	char *rhofile="";	/* density filename*/
	FILE *rhofilefp=NULL;	/* ... its file pointer */
	char *c11file="";	/* c11 stiffness filename */
	FILE *c11filefp=NULL;	/* ... its file pointer */
	char *c55file="";	/* c55 stiffness filename */
	FILE *c55filefp=NULL;	/* ... its file pointer */
	char *c13file="";	/* c13 stiffness filename */
	FILE *c13filefp=NULL;	/* ... its file pointer */
	char *c15file="";	/* c15 stiffness filename */
	FILE *c15filefp=NULL;	/* ... its file pointer */
	char *c33file="";	/* c33 stiffness filename */
	FILE *c33filefp=NULL;	/* ... its file pointer */
	char *c35file="";	/* c35 stiffness filename */
	FILE *c35filefp=NULL;	/* ... its file pointer */
	char *qfile="";	 /* Q filename*/
	FILE *qfilefp=NULL;     /* ... its file pointer */


	float sa,ca,p;

	int mode;		/* output mode =0 for particle velocity */
				/* =1 for stresses */

	/* hook up getpar to handle the parameters */
	initargs(argc,argv);
	requestdoc(0);

	/* get optional parameters */
	if (!getparfloat("dt",&dt)) dt = 0.001;
	if (!getparfloat("ft",&ft)) ft = 0.0;
	if (!getparfloat("lt",&lt)) lt = 1.0;

	/* get model size in slow direction */
	if (!getparint("nx",&nx)) nx = 200;
	if (!getparfloat("dx",&dx)) dx = 10;
	if (!getparfloat("fx",&fx)) fx = -1000;

	/* get model size in fast direction */
	if (!getparint("nz",&nz)) nz = 100;
	if (!getparfloat("dz",&dz)) dz = dx;
	if (!getparfloat("fz",&fz)) fz = 0;

	/* calculate maximum model coordinates */
	xmax = fx + dx*nx;
	zmax = fz + dz*nz;

	/* get source location */
	if (!getparfloat("sx",&sx)) sx = (fx+xmax)/2;
	if (!getparfloat("sz",&sz)) sz = (fz+zmax)/2;

	/* get source type, angle */	
	if (!getparstring("stype",&stype)) stype="p";
	if (!getparfloat("sang",&sang)) sang = 0;

	/* get wavelet type, source time length, average frequency */
	if (!getparstring("wtype",&wtype)) wtype="ri";
	if (!getparfloat("ts",&ts)) ts = 0.05;
	if (!getparfloat("favg",&favg)) favg = 50;

	if (!getparint("qsw",&qsw)) qsw = 0;
	if (!getparint("aniso",&aniso)) aniso = 0;
	if (!getparint("tsw",&tsw)) tsw = 0;
	if (!getparint("verbose",&verbose)) verbose = 0;
	if (!getparint("mode",&mode)) mode = 0;

	/***...begin......... shapshot and source files....... */
	/* get name of output snapshot, source files */ 
	getparstring("snfile",&snfile);
	getparstring("sofile",&sofile);

	/* determine snapshots times and sample numbers */
	nsnap = countparval("snaptime");
	if (nsnap!=0){
		snaptime = alloc1float(nsnap);
		isnap = alloc1int(nsnap);
		getparfloat("snaptime",snaptime);
		warn("nsnap=%d \n",nsnap);
		for (i=0; i<nsnap; i++) isnap[i]=NINT(snaptime[i]/dt);
	}
	/* if requested, open file for snapshots */
	if (*snfile!='\0') {
		if((sneisfp=fopen(snfile,"w"))==NULL)
			err("cannot open snfile=%s\n",snfile);
	}

	/* if requested, open file for source */
	if (*sofile!='\0') {
		if((soeisfp=fopen(sofile,"w"))==NULL)
			err("cannot open sofile=%s\n",sofile);
	}
	/***...end......... shapshot and source files....... */


	/**** begin ....... stiffnesses  .... */
	/* Must have stiffnesses c11 and c55 */
	if (!getparstring("c11file",&c11file)) 	c11file="c11_file";
	if (!getparstring("c55file",&c55file)) 	c55file="c55_file";

	/*.. open file for C11 */
	if((c11filefp=fopen(c11file,"r"))==NULL)
			err("cannot open c11file=%s\n",c11file);
	/*.. open file for C55 */
	if((c55filefp=fopen(c55file,"r"))==NULL)
			err("cannot open c55file=%s\n",c55file);


	if (aniso==1) { 
		/* Must have these additional stiffnesses for anisotropy */
		if (!getparstring("c13file",&c13file))	c13file="c13_file";
		if (!getparstring("c35file",&c35file))	c35file="c35_file";
		if (!getparstring("c33file",&c33file))	c33file="c33_file";
		if (!getparstring("c15file",&c15file))	c15file="c15_file";

		/*.. open file for C13 */
		if((c13filefp=fopen(c13file,"r"))==NULL)
				err("cannot open c13file=%s\n",c13file);

		/*.. open file for C35 */
		if((c35filefp=fopen(c35file,"r"))==NULL)
				err("cannot open c35file=%s\n",c35file);

		/*.. open file for C33 */
		if((c33filefp=fopen(c33file,"r"))==NULL)
				err("cannot open c33file=%s\n",c33file);

		/*.. open file for C15 */
		if((c15filefp=fopen(c15file,"r"))==NULL)
				err("cannot open c15file=%s\n",c15file);
	}
	/**** end ....... stiffnesses .... */

	/**** start ..... attenuation .... */
	if (qsw==1) {
		/* Must have these values for attenuation */
		if (!getparstring("qfile",&qfile))	qfile="q_file";

		/*.. open file for Q */
		if((qfilefp=fopen(qfile,"r"))==NULL)
				err("cannot open qfile=%s\n",qfile);
	}
	/**** end ....... attenuation .... */

	
	/**** begin ..... boundary conditions .... */
	/* get boundary conditions */
	nbc = countparval("bc");
	if (nbc==4) {
		getparint("bc", bc);	
	} else {
		bc[0] = 10;
		bc[1] = 10;
		bc[2] = 10;
		bc[3] = 10;

		if (!((nbc==4) || (nbc==0)) ) 
			warn("Number of bc %d, using bc=10,10,10,10",bc);
	}
	if (!getparfloat("bc_a",&bc_a)) bc_a = 0.95;
	if (!getparfloat("bc_r",&bc_r)) bc_r = 0;
	warn("%f bc_a  %f bc_r\n",bc_a,bc_r);

	for(i=0; i<4; i++) {
		wbc[i]=2;
		if(bc[i]>2)wbc[i]=bc[i];
	}

	bc0=alloc1float(bc[0]);
	bc0[0]=bc_a;
	bc1=alloc1float(bc[1]);
	bc1[0]=bc_a;
	bc2=alloc1float(bc[2]);
	bc2[0]=bc_a;
	bc3=alloc1float(bc[3]);
	bc3[0]=bc_a;

	for(j=1; j<bc[0]; j++) bc0[j]=bc0[j-1]*bc_a*pow(j,-bc_r);
	for(j=1; j<bc[1]; j++) bc1[j]=bc1[j-1]*bc_a*pow(j,-bc_r);
	for(j=1; j<bc[2]; j++) bc2[j]=bc2[j-1]*bc_a*pow(j,-bc_r);
	for(j=1; j<bc[3]; j++) bc3[j]=bc3[j-1]*bc_a*pow(j,-bc_r);

		/*warn("%d %f\n",j,bc0[j]);*/

 	for(j=0; j<bc[3]; j++) warn(" %d %f\n",j,bc3[j]);
	
	warn("bc=%d %d %d %d",bc[0],bc[1],bc[2],bc[3]);
	warn("wbc=%d %d %d %d",wbc[0],wbc[1],wbc[2],wbc[3]);
	/**** end ..... boundary conditions .... */
	
	/***** begin .... seismogram files .... */
	/* determine seismogram recording lines */
	if (!getparfloat("hsz",&hsz)) hsz = 0.0;
	hs1 = NINT( (hsz - fz)/dz )+wbc[0];		/*horizontal line*/
	if (!getparfloat("vsx",&vsx)) vsx = 0.0;
	vs2 = NINT((vsx - fx)/dx )+wbc[1];		/*vertical line*/
	if (!getparint("verbose",&verbose)) verbose = 0;
	if (!getparstring("hsfile",&hsfile)) hsfile="hs.su";
	if (!getparstring("vrslfile",&vrslfile)) vrslfile="vsp.su";
	warn("hs1=%d vs2=%d \n",hs1, vs2);
	

	/* open files for seismograms */
	if((hseisfp=fopen(hsfile,"w+"))==NULL) 
		err("cannot open hsfile=%s\n",hsfile);
	if((vseisfp=fopen(vrslfile,"w+"))==NULL)
		err("cannot open vrslfile=%s\n",vrslfile);
	/***** end .... seismogram files .... */
		
	/* pad the nx and nz dimensions for boundary condition layers */
	nxpadded=nx+wbc[1]+wbc[3];		
	nzpadded=nz+wbc[0]+wbc[2];

	nt=NINT((lt-ft)/dt)+1;
	warn("dx=%f dz=%f dt=%f",dx,dz,dt);
	warn("nxpadded=%d nzpadded=%d nt=%d",nxpadded,nzpadded,nt);
	warn("array sizes");

	/* set grid parameters*/
	limits = alloc1int(4);
	 	dtx=dt/dx; dtz=dt/dz;
	 	energy=0;efac=1000000000000000000./((nxpadded-4)*(nzpadded-4));
	 	pfac=1./(dx*dz);
	 	is=NINT((-fx+sx)/dx)+wbc[1]; js=NINT((-fz+sz)/dz)+wbc[0];
 	warn("is=%d js=%d\n",is,js);

	/* check to se if source is inside the model */
 	if((sx<fx) || (sx>xmax))
		err("sx=%f source must lie within model \n",sx);
 	if((sz<fz) || (sz>zmax))
		err("sz=%f source must lie within model \n",sz);
		
	/* allocate space for model*/
	c11 = alloc2float(nxpadded,nzpadded);
	c11temp = alloc2float(nx,nz);
	c55 = alloc2float(nxpadded,nzpadded);
	c55temp = alloc2float(nx,nz);
	rho = alloc2float(nxpadded,nzpadded);
	rhotemp = alloc2float(nx,nz);
	ngrids=3;

	/* if attenuation is desired allocate space for q values */
	if(qsw==1) {
		q = alloc2float(nxpadded,nzpadded);
		qtemp = alloc2float(nx,nz);
		++ngrids; /* bump ngrids up by one */
	}

	/* calculate the number of bytes in grids */
	nbytes=ngrids*nxpadded*nzpadded*FSIZE;

	/* if anisotropy is desired allocate space for anisotropy params */
	if(aniso==1) {
		c33 = alloc2float(nxpadded,nzpadded);
		c33temp = alloc2float(nx,nz);
		c13 = alloc2float(nxpadded,nzpadded);
		c13temp = alloc2float(nx,nz);
		c15 = alloc2float(nxpadded,nzpadded);
		c15temp = alloc2float(nx,nz);
		c35 = alloc2float(nxpadded,nzpadded);
		c35temp = alloc2float(nx,nz);
		ngrids+=4;
	}

	/* calculate the number of bytes in all the grids */
	nbytes=ngrids*nxpadded*nzpadded*FSIZE;
	warn("%d grids allocated for model",ngrids);
	warn("%f MBytes allocated for model",0.000001*nbytes);

	/* set e_*/
	strcpy(tfile,hsfile);
	mfile=strcat(tfile,".mod");

	/* open mfile */	
	if((mfp=fopen(mfile,"w+"))==NULL)
		err("cannot open mfile=%s\n",mfile);

	/***** get density and stiffnesses from input files ***********/
	/* ... densities optional */
	getparstring("rhofile",&rhofile);

 	/*** begin density **/
	/* if specified, read density rho */
	if (*rhofile!='\0') {
		if((rhofilefp=fopen(rhofile,"r"))==NULL)
			err("cannot open rhofile=%s\n",rhofile);
		if (fread(rhotemp[0],sizeof(float),nx*nz,rhofilefp)!=nx*nz)
			err("error reading rhofile=%s\n",rhofile);
		fclose(rhofilefp);

		/* pad the density */
		pad_float_2D(nx,nz,wbc,rhotemp,rho); 
		free2float(rhotemp);



		/* we need to invert density */
		/* find min and max of 1/density */
		rhomin = rhomax = 1/rho[0][0];
		for (iz=0; iz<nzpadded; ++iz) {
			for (ix=0; ix<nxpadded; ++ix) {

					/* rho = 1/rho */
					rho[iz][ix] = 1/rho[iz][ix];

					rhomin= MIN(rhomin,rho[iz][ix]);
					rhomax = MAX(rhomax,rho[iz][ix]);
			}
		}
			
	 }

	/* if densities not specified  make densities = 10^{-3} */
	if (*rhofile=='\0' ) {
		err("!!!rhofile= not set!!!");
	  /* if (verbose) warn("!!!rhofile= not set, assuming density=1000!!!");
		  for (iz=0; iz<nzpadded; ++iz)
			   for (ix=0; ix<nxpadded; ++ix)
				    rho[iz][ix] = .001;
		  rhomin = rhomax = .001 ;
	  */
	}
	if (verbose) warn(" 1/density:  max = %f  min=%f",rhomax,rhomin);
 	/** end density **/


	/***** begin c11  */
	if((c11filefp=fopen(c11file,"r"))==NULL)
		err("cannot open c11file=%s\n",c11file);
	if (fread(c11temp[0],sizeof(float),nx*nz,c11filefp)!=nx*nz)
		err("error reading c11file=%s\n",c11file);
	fclose(c11filefp);

	/* pad the c11 */
	pad_float_2D(nx,nz,wbc,c11temp,c11); 
	free2float(c11temp);

	/* find min and max of c11 */
	c11min = c11max = c11[0][0];
	for (iz=0; iz<nzpadded; ++iz) {
		for (ix=0; ix<nxpadded; ++ix) {
			    c11min = MIN(c11min,c11[iz][ix]);
			    c11max = MAX(c11max,c11[iz][ix]);
		}
	}
	if (verbose)  warn(" c11max=%f  c11min=%f",c11max,c11min);
			
 	/**** end ....... c11 **/

	/***** begin ...........c55  */
	if((c55filefp=fopen(c55file,"r"))==NULL)
		err("cannot open c55file=%s\n",c55file);
	if (fread(c55temp[0],sizeof(float),nx*nz,c55filefp)!=nx*nz)
		err("error reading c55file=%s\n",c55file);
	fclose(c55filefp);

	/* pad the c55 */
	pad_float_2D(nx,nz,wbc,c55temp,c55); 
	free2float(c55temp);

	/* find min and max of c55 */
	c55min = c55max = c55[0][0];
	for (iz=0; iz<nzpadded; ++iz) {
		for (ix=0; ix<nxpadded; ++ix) {
			    c55min = MIN(c55min,c55[iz][ix]);
			    c55max = MAX(c55max,c55[iz][ix]);
		}
	}
	if (verbose) warn(" c55max=%f  c55min=%f",c55max,c55min);

			
 	/**** end ....... c55 **/

	if (aniso==1) {  /* if anisotropy */

		/***** begin ...........c13  */
		if((c13filefp=fopen(c13file,"r"))==NULL)
			err("cannot open c13file=%s\n",c13file);
		if (fread(c13temp[0],sizeof(float),nx*nz,c13filefp)!=nx*nz)
			err("error reading c13file=%s\n",c13file);
		fclose(c13filefp);

		/* pad the c13 */
		pad_float_2D(nx,nz,wbc,c13temp,c13); 
		free2float(c13temp);

		/* find min and max of c13 */
		c13min = c13max = c13[0][0];
		for (iz=0; iz<nzpadded; ++iz) {
			for (ix=0; ix<nxpadded; ++ix) {
				    c13min = MIN(c13min,c13[iz][ix]);
				    c13max = MAX(c13max,c13[iz][ix]);
			}
		}
			
 		/**** end ....... c13 **/

		/***** begin ...........c33  */
		if((c33filefp=fopen(c33file,"r"))==NULL)
			err("cannot open c33file=%s\n",c33file);
		if (fread(c33temp[0],sizeof(float),nx*nz,c33filefp)!=nx*nz)
			err("error reading c33file=%s\n",c33file);
		fclose(c33filefp);

		/* pad the c33 */
		pad_float_2D(nx,nz,wbc,c33temp,c33); 
		free2float(c33temp);

		/* find min and max of c33 */
		c33min = c33max = c33[0][0];
		for (iz=0; iz<nzpadded; ++iz) {
			for (ix=0; ix<nxpadded; ++ix) {
				    c33min = MIN(c33min,c33[iz][ix]);
				    c33max = MAX(c33max,c33[iz][ix]);
			}
		}
			
 		/**** end ....... c33 **/

		/***** begin ...........c15  */
		if((c15filefp=fopen(c15file,"r"))==NULL)
			err("cannot open c15file=%s\n",c15file);
		if (fread(c15temp[0],sizeof(float),nx*nz,c15filefp)!=nx*nz)
			err("error reading c15file=%s\n",c15file);
		fclose(c15filefp);

		/* pad the c15 */
		pad_float_2D(nx,nz,wbc,c15temp,c15); 
		free2float(c15temp);

		/* find min and max of c15 */
		c15min = c15max = c15[0][0];
		for (iz=0; iz<nzpadded; ++iz) {
			for (ix=0; ix<nxpadded; ++ix) {
				    c15min = MIN(c15min,c15[iz][ix]);
				    c15max = MAX(c15max,c15[iz][ix]);
			}
		}
 		/**** end ....... c15 **/

		/***** begin ...........c35  */
		if((c35filefp=fopen(c35file,"r"))==NULL)
			err("cannot open c35file=%s\n",c35file);
		if (fread(c35temp[0],sizeof(float),nx*nz,c35filefp)!=nx*nz)
			err("error reading c35file=%s\n",c35file);
		fclose(c35filefp);

		/* pad the c35 */
		pad_float_2D(nx,nz,wbc,c35temp,c35); 
		free2float(c35temp);

		/* find min and max of c35 */
		c35min = c35max = c35[0][0];
		for (iz=0; iz<nzpadded; ++iz) {
			for (ix=0; ix<nxpadded; ++ix) {
				    c35min = MIN(c35min,c35[iz][ix]);
				    c35max = MAX(c35max,c35[iz][ix]);
			}
		}
			
 		/**** end ....... c35 **/
	}

	if (qsw==1) {  /* if attenuation */

		if((qfilefp=fopen(qfile,"r"))==NULL)
			err("cannot open qfile=%s\n",qfile);
		if (fread(qtemp[0],sizeof(float),nx*nz,qfilefp)!=nx*nz)
			err("error reading qfile=%s\n",qfile);
		fclose(qfilefp);

		/* pad q */
		pad_float_2D(nx,nz,wbc,qtemp,q);
		free2float(qtemp);

		/* find min and max of q */
		qmin = qmax = q[0][0];
		for (iz=0; iz<nzpadded; ++iz) {
			for (ix=0; ix<nxpadded; ++ix) {
				    qmin = MIN(qmin,q[iz][ix]);
				    qmax = MAX(qmax,q[iz][ix]);
			}
		}
		if (verbose) warn(" qmax=%f  qmin=%f",qmax,qmin);
	}


	if (verbose) warn(" c11max=%f  c11min=%f",c11max,c11min);
	if (verbose) warn(" rhomax=%f  rhomin=%f",rhomax,rhomin);
	if (verbose) warn(" wbc[0]=%d  wbc[1]=%d wbc[2]=%d wbc[3]=%d",
			wbc[0],wbc[1],wbc[2],wbc[3]);

	vlim[0] = sqrt(c11min*rhomax);
	vlim[1] = sqrt(c11max*rhomin);


	/* close mfile */
	fclose(mfp);

	/* set arrays for plane boundary conditions */
	ttl = alloc1float(nzpadded); 
	ttb = alloc1float(nxpadded);
	ttr = alloc1float(nzpadded);
	vl = alloc1float(nzpadded);
	vb = alloc1float(nxpadded);
	vr = alloc1float(nzpadded);

	/* note that rho = 1/density now */
	 for (j=0; j<nzpadded; ++j) {
		vl[j]=sqrt(c11[j][wbc[1]]*rho[j][wbc[1]]);
		vr[j]=sqrt(c11[j][nxpadded-wbc[3]-1]*rho[j][nxpadded-wbc[3]-1]);
	}
	 	for (i=0; i<nxpadded; ++i) {
		vb[i]=sqrt(c11[nzpadded-wbc[2]-1][i]*rho[nzpadded-wbc[2]-1][i]);
	}

	/* ttb: traveltime from source to location on bottom */
	sa=sin(PI*sang/180.);ca=cos(PI*sang/180.);

	if(sang>=0&&sang<90){
		ttb[0]=(-wbc[1])*sa*dx/vb[wbc[1]];
	 		for (i=1; i<nxpadded; ++i) {
			ttb[i]=ttb[i-1]+dx*sa/vb[i];
		}
	}

	if(sang>-90&&sang<0){
		ttb[nxpadded-1]=-(-wbc[3])*sa*dx/vl[nxpadded-wbc[3]-1];
	 		for (i=nxpadded-1; i>0; --i) {
			ttb[i-1]=ttb[i]-dx*sa/vb[i];
		}
	}

	/* ttl, ttr: traveltimes from source to location on LHS, RHS */
	sa=sin(PI*sang/180.);
	ca=cos(PI*sang/180.);
	p=sin(sang*PI/180.)/vl[nzpadded-wbc[2]-1];

	ttl[nzpadded-1]=(-wbc[2])*ca*dz/vl[nzpadded-wbc[2]-1];
	 	for (j=nzpadded-1; j>0; --j) {
		sa=sin(asin(p*vl[j]));
		ca=cos(asin(p*vl[j]));
		ttl[j-1]=ttl[j]+dz*ca/vl[j];
	}

	sa=sin(PI*sang/180.);
	ca=cos(PI*sang/180.);
	p=sin(sang*PI/180.)/vr[nzpadded-wbc[2]-1];
	ttr[nzpadded-1]=(-wbc[2])*ca*dz/vr[nzpadded-wbc[2]-1];

	 	for (j=nzpadded-1; j>0; --j) {
		sa=sin(asin(p*vr[j]));
		ca=cos(asin(p*vr[j]));
		ttr[j-1]=ttr[j]+dz*ca/vr[j];
	}

	warn("vmin=%f vmax=%f",vlim[0],vlim[1]);
	warn("model set");
	
	/* calculate stability and dispersion */
	warn("for stablity dt should be < %f: x-dir and %f: z-dir",
		0.6*dx/vlim[1],0.6*dz/vlim[1]);
	warn("for low dispersion fmax should be < %f ",
		vlim[0]/(5*dx));
	
		
	/* get source function */
	ns=NINT(ts/dt)+1;
	source = alloc1float(ns);
	i=get_source(dt,ts,favg,wtype,source);
	warn("source set");
	if (*sofile!='\0') {
		fprintf(soeisfp,"%s %s %d %d %f %f\n",wtype,stype,ns,i,favg,ts);
		for (i=0; i<ns; i++)
			fprintf(soeisfp,"%f %f\n",i*dt,source[i]);
 	}

	pfac=1.;
	for (i=0; i<ns; i++) source[i]*=pfac;

 	warn("%s %s %d %d %f %f",wtype,stype,ns,i,favg,ts);
	
	/* allocate space for variables*/
	/*  CJ: ADJUSTMENT MADE 030502 - txx,tzz and w extended for symmetry */
	nbytes=ngrids*nxpadded*nzpadded*FSIZE;
	u = alloc2float(nxpadded,nzpadded);
	nbytes=nbytes+nxpadded*nzpadded*FSIZE;
	w = alloc2float(nxpadded+1,nzpadded);
	nbytes=nbytes+(nxpadded+1)*nzpadded*FSIZE;
	txx = alloc2float(nxpadded+1,nzpadded);
	nbytes=nbytes+(nxpadded+1)*nzpadded*FSIZE;
	tzz = alloc2float(nxpadded+1,nzpadded);
	nbytes=nbytes+(nxpadded+1)*nzpadded*FSIZE;
	txz = alloc2float(nxpadded,nzpadded);
	nbytes=nbytes+nxpadded*nzpadded*FSIZE;
	ngrids+=5;
	
	warn("%f MBytes allocated for model + variables\n",
			0.000001*nbytes);

	/* scale Q-factor */
	if(qsw==1) {
		for (i=0; i<nxpadded; i++){
			for (j=0; j<nzpadded; j++) {
				q[j][i]=exp(-PI*dt*favg/q[j][i]); 
			}
		}	
	}
	
	/* zero arrays */
	for (i=0; i<nxpadded; i++){
		for (j=0; j<nzpadded; j++) {
			u[j][i]=0;
			w[j][i]=0;
			txx[j][i]=0;
			tzz[j][i]=0;
			txz[j][i]=0;
		}
	}

	/* write character headers */
	strcpy(tfile,hsfile);
	hcfile=strcat(tfile,".chd");
	if((hchdfp=fopen(hcfile,"w+"))==NULL)
		err("cannot open hcfile=%s\n",hcfile);
	write_chd(nxpadded,nzpadded,nt,dx,dz,dt,fx,xmax,fz,zmax,
			sx,sz,favg,ts,wtype,stype,
			bc,qsw,aniso,hsz,hchdfp,1,mfile);
	fclose(hchdfp);

	strcpy(tfile,vrslfile);vrshfile=strcat(tfile,".chd");
	if((vchdfp=fopen(vrshfile,"w+"))==NULL)
		err("cannot open vrshfile=%s\n",vrshfile);
	write_chd(nxpadded,nzpadded,nt,dx,dz,dt,fx,xmax,fz,zmax,
			sx,sz,favg,ts,wtype,stype,
			bc,qsw,aniso,vsx,vchdfp,2,mfile);
	fclose(vchdfp);

	/* begin wave propagation */
	k=0;
	for (t=ft; t<=lt; t=t+dt) {

		/* calculate area for active grid */
		calc_area(vlim[1],dtx,dtz,limits,is,js,k,nxpadded,nzpadded);

		ifx=limits[0];
		ilx=limits[1];
		jfz=limits[2];
		jlz=limits[3];

		if (strcmp(stype,"pw")==0){
			ifx=2;
			ilx=nxpadded-3;
			jfz=2;
			jlz=nzpadded-3;
		}

		/* introduce velocity source */
		if (k<ns) {
			if (strcmp(stype,"v")==0)
				add_v_source(u,w,dt*source[k],is,js);
		}

		/* update velocities */
		if(tsw==0)
			update_vel(ifx,ilx,jfz,jlz,u,w,txx,tzz,txz,
					rho,c55,dtx,dtz);
		if(tsw==1)
			update_vel_tsw(ifx,ilx,jfz,jlz,u,w,txx,tzz,
					txz,rho,c55,dtx,dtz);
	
		/* introduce plane wave velocity source */
		if (strcmp(stype,"pw")==0){
		  		add_pw_source_V(u,w,txx,tzz,txz,c11,c55,
						source,sang,wbc[1],
						nxpadded-wbc[3]-1,
						wbc[0],nzpadded-wbc[2]-1,
						ns,k,dx,dz,dt,vl,vb,vr,
						ttl,ttb,ttr);
		}

		/* apply free surface bc for velocities */
		if(bc[0]==2) fs4v_bc_top(u,w,c11,c55,nxpadded,nzpadded);
	
		/* apply symmetry bc for velocities */
		if(bc[0]==1) sym4v_bc_top(u,w,nxpadded,nzpadded);
		if(bc[1]==1) sym4v_bc_left(u,w,nxpadded,nzpadded);
		if(bc[2]==1) sym4v_bc_bot(u,w,nxpadded,nzpadded);
		if(bc[3]==1) sym4v_bc_right(u,w,nxpadded,nzpadded);
	
		/* introduce stress source */
		if (k<ns) {
	  		if (strcmp(stype,"p")==0)
				add_p_source(txx,tzz,dt*source[k],is,js);
		}
	
		/* update stresses */
		if(!aniso) {  /* isotropic */
			update_stress_iso(ifx,ilx,jfz,jlz,u,w,txx,tzz,txz,
						c11,c55,dtx,dtz);
		} else {
			update_stress_ani(ifx,ilx,jfz,jlz,u,w,txx,tzz,txz,
						c11,c55,c33,c13,c15,c35,
						dtx,dtz);
		}

		/* introduce plane wave stress source */
		if (strcmp(stype,"pw")==0) {
			add_pw_source_S(u,w,txx,tzz,txz,c11,c55,
						source,sang,wbc[1],
						nxpadded-wbc[3]-1,wbc[0],
						nzpadded-wbc[2]-1,ns,k,
						dx,dz,dt,vl,vb,vr,ttl,
						ttb,ttr);
		}

		/* apply free surface bc for stresses */
		if(bc[0]==2) fs4s_bc_top(txx,tzz,txz,nxpadded,nzpadded);

		/* apply symmetry bc for stresses */
		if(bc[0]==1) sym4s_bc_top(txx,tzz,txz,nxpadded,nzpadded);
		if(bc[1]==1) sym4s_bc_left(txx,tzz,txz,nxpadded,nzpadded);
		if(bc[2]==1) sym4s_bc_bot(txx,tzz,txz,nxpadded,nzpadded);
		if(bc[3]==1) sym4s_bc_right(txx,tzz,txz,nxpadded,nzpadded);

		/* apply absorbing bc */
		if(bc[0]>2)
			abs_bc_top(u,w,txx,tzz,txz,nxpadded,nzpadded,
					bc0,bc,wbc,ifx,ilx);
		if(bc[1]>2)
			abs_bc_left(u,w,txx,tzz,txz,nxpadded,nzpadded,
					bc1,bc,wbc,jfz,jlz);
		if(bc[2]>2)
			abs_bc_bot(u,w,txx,tzz,txz,nxpadded,nzpadded,
					bc2,bc,wbc,ifx,ilx);
		if(bc[3]>2)
			abs_bc_right(u,w,txx,tzz,txz,nxpadded,nzpadded,
					bc3,bc,wbc,jfz,jlz);

		/* attenuation */
		if(qsw==1) {
			for (j=jfz-2; j<jlz+2; j++){
				for (i=ifx-2; i<ilx+2; i++) {
			 		u[j][i]*=q[j][i];
					w[j][i]*=q[j][i];
			 		txx[j][i]*=q[j][i];
					tzz[j][i]*=q[j][i];
					txz[j][i]*=q[j][i];
				}
			}	
		}
	
		/* "energy" */
		energy=0;
		for (i=ifx; i<ilx; i++){
			for (j=jfz; j<jlz; j++) {
				energy=energy+u[j][i]*u[j][i]+w[j][i]*w[j][i]; 
			}
		}
	
		/* store wave motion on file (seismograms) */
		for (i=0 ; i < nxpadded ; ++i)
			fwrite(&u[hs1][i], 4, 1, hseisfp);
		for (i=0 ; i < nxpadded ; ++i) 
			fwrite(&w[hs1][i], 4, 1, hseisfp);
		for (j=0 ; j < nzpadded ; ++j) 
			fwrite(&u[j][vs2], 4, 1, vseisfp);
		for (j=0 ; j < nzpadded ; ++j) 
			fwrite(&w[j][vs2], 4, 1, vseisfp);
				
		/* make snapshot */
		for (i=0; i<nsnap; i++) {
			if (!mode) {
				if(k==isnap[i] && *snfile!='\0') {
					make_snap(u,w,sx,sz,dx,dz,
							nxpadded,nzpadded,
							nx,nz,
							k,t,fx,fz,tr,
							sneisfp,wbc);
		 		} 
			} else {
				if(k==isnap[i] && *snfile!='\0') {
		  			make_stress_snap(txx,tzz,txz,sx, sz,
								dx, dz,
								nxpadded, 
								nzpadded,
								k,t,fx,fz,tr,
								sneisfp,wbc);
				}
		 	}
	 	}

		if((verbose==1) && (k%50==0))
			fprintf(stderr,"%d %f %f %d %d %d %d \n",
 				k,t,energy*efac,ifx,ilx,jfz,jlz);


		k=k+1;
	}

	warn("propagation completed\n");

	/* free space */
	free2float(u);
	free2float(w);
	free2float(txx);
	free2float(tzz);
	free2float(txz);
	warn("propagation space freed\n");
	
	/* read in wave motion from file and make seismograms */
	rewind(hseisfp);	/*horizontal line*/ 

	ut = alloc2float(nxpadded,nt); 
	wt = alloc2float(nxpadded,nt);

	for (k=0 ; k < nt ; ++k){
		for (i=0 ; i < nxpadded ; ++i) 
			fread(&ut[k][i], 4, 1, hseisfp);
		for (i=0 ; i < nxpadded ; ++i)
			fread(&wt[k][i], 4, 1, hseisfp);
	}
	fclose(hseisfp);
	hseisfp=fopen(hsfile,"w");

	make_seis(ut,wt,sx,sz,dx,dt,nxpadded,nt,hsz,fx,fz,
			trh,hseisfp,wbc[1],wbc[3],1,verbose,favg);
	fclose(hseisfp);
	free2float(ut);
	free2float(wt);
	warn("u time section output\n");
	
	rewind(vseisfp);	/*vertical line*/
	ut = alloc2float(nzpadded,nt); wt = alloc2float(nzpadded,nt);
	for (k=0 ; k < nt ; ++k){
	  for (j=0 ; j < nzpadded ; ++j) fread(&ut[k][j], 4, 1, vseisfp);
	  for (j=0 ; j < nzpadded ; ++j) fread(&wt[k][j], 4, 1, vseisfp);
	}
	fclose(vseisfp); vseisfp=fopen(vrslfile,"w");

	make_seis(ut,wt,sx,sz,dz,dt,nzpadded,nt,vsx,fx,fz,
			trv,vseisfp,wbc[0],wbc[2],2,verbose,favg);

	fclose(vseisfp); free2float(ut); free2float(wt);
	warn("w time section output\n");

	if (mode!=0) fclose(sneisfp);

	return(CWP_Exit());
}

int get_source(float dt, float ts, float favg, char *wtype, float *source)
{
	int i,ns; float t,x,xx,pi2;
	i=0; t=0; pi2=PI*PI;
	ns=NINT(ts/dt)+1;
	while (i<ns) {
	  if(strcmp(wtype,"ri")==0){ /*Ricker*/
		x=favg*(t-ts/2);xx=x*x;
		source[i]=(1-2*pi2*(xx))*exp(-pi2*xx);
	  } 							  
	  if(strcmp(wtype,"dg")==0){ /*Derivative of Gaussian*/
		x=-4*favg*favg*pi2/log(0.1);
		source[i]=-2*pi2*(t-ts/2)*exp(-x*(t-ts/2)*(t-ts/2));
	  } 							  
	  if(strcmp(wtype,"ga")==0){ /*Gaussian favg is 10% level for fmax*/
		x=-favg*favg*pi2/log(0.1);
		source[i]=exp(-x*(t-ts/2)*(t-ts/2));
	  } 							  
	  if(strcmp(wtype,"sp")==0){ /*spike*/
		source[0]=1;
	  } 						
	  if(strcmp(wtype,"sp2")==0){ /*spike*/
		source[0]=1;
		source[1]=1;
	  } 						
 	  /*warn("%d %f %f\n",i,t,source[i]);*/
	  t=t+dt; i=i+1;
	}
	return(i);
}

void update_vel(int ifx, int ilx, int jfz, int jlz,
	float **u, float **w, float **txx, float **tzz, float **txz,
	float **rho, float **c55, float dtx, float dtz)
{
	int i,j;
	float dtxx=0,dtzz=0,dtxz=0,dtzx=0;
	 	for (j=jfz; j<=jlz; ++j) {
	 	  for (i=ifx; i<=ilx; ++i) {
	 		dtxx = F1*(txx[j][i+1]-txx[j][i])
		 +F2*(txx[j][i+2]-txx[j][i-1]);
		dtxz = F1*(txz[j+1][i]-txz[j][i])
	 		+F2*(txz[j+2][i]-txz[j-1][i]);
			u[j][i] = u[j][i] + ( dtx*dtxx + dtz*dtxz ) * rho[j][i];
	 	  }
	 	  for (i=ifx; i<=ilx+1; ++i) {
			dtzz = F1*(tzz[j][i]-tzz[j-1][i])
	 		+F2*(tzz[j+1][i]-tzz[j-2][i]);
	 		dtzx = F1*(txz[j][i]-txz[j][i-1])
	 		+F2*(txz[j][i+1]-txz[j][i-2]);
	 		w[j][i] = w[j][i] + ( dtx*dtzx + dtz*dtzz ) * rho[j][i];
	 	  }
	 	}
}
void update_vel_tsw(int ifx, int ilx, int jfz, int jlz,
	float **u, float **w, float **txx, float **tzz, float **txz,
	float **rho, float **c55, float dtx, float dtz)
{
	int i,j;
	float dtxx=0,dtzz=0,dtxz=0,dtzx=0;
	 	for (j=jfz; j<=jlz; ++j) {
	 	  for (i=ifx; i<=ilx; ++i) {
	 		dtxx = F1*(txx[j][i+1]-txx[j][i])
		 +F2*(txx[j][i+2]-txx[j][i-1]);
/*  		Only use 4th order shear stress derivatives if material is non-fluid, otherwise problems occur at grazing angles near fluid-solid boundary */
	 		if(c55[j][i]>0) { 
		 dtxz = F1*(txz[j+1][i]-txz[j][i])
	 		+F2*(txz[j+2][i]-txz[j-1][i]);
		}
		else{
		 dtxz = (txz[j+1][i]-txz[j][i]);
		}
	 		u[j][i] = u[j][i] + ( dtx*dtxx + dtz*dtxz ) * rho[j][i];
	 	  }
	 	  for (i=ifx; i<=ilx+1; ++i) {
	 		dtzz = F1*(tzz[j][i]-tzz[j-1][i])
	 		+F2*(tzz[j+1][i]-tzz[j-2][i]);
/*  		Only use 4th order shear stress derivatives if material is non-fluid, otherwise problems occur at grazing angles near fluid-solid boundary */
	 		if(c55[j][i]>0) { 
	 		dtzx = F1*(txz[j][i]-txz[j][i-1])
	 		+F2*(txz[j][i+1]-txz[j][i-2]);
		}
		else{
	 		dtzx = (txz[j][i]-txz[j][i-1]);
		}
	 		w[j][i] = w[j][i] + ( dtx*dtzx + dtz*dtzz ) * rho[j][i];
	 	  }
	 	}
}


void update_stress_iso(int ifx, int ilx, int jfz, int jlz,
	float **u, float **w, float **txx, float **tzz, float **txz,
	float **c11, float **c55, float dtx, float dtz)
{
	int i,j;
	float dux,duz,dwz,dwx;
	 	for (j=jfz; j<=jlz; ++j) {
	 	  for (i=ifx; i<=ilx+1; ++i) {
	 		dwz = F1*(w[j+1][i]-w[j][i])
	 		+F2*(w[j+2][i]-w[j-1][i]);
	 		dux = F1*(u[j][i]-u[j][i-1])
	 		+F2*(u[j][i+1]-u[j][i-2]);
	 		txx[j][i] = txx[j][i] +
	 		c11[j][i]*dtx*dux + (c11[j][i]-2*c55[j][i])*dtz*dwz;
	 		tzz[j][i] = tzz[j][i] +
	 		c11[j][i]*dtz*dwz + (c11[j][i]-2*c55[j][i])*dtx*dux;
	 	  }
	 	  for (i=ifx; i<=ilx; ++i) {
	 		dwx = F1*(w[j][i+1]-w[j][i])
		 +F2*(w[j][i+2]-w[j][i-1]);
	 		duz = F1*(u[j][i]-u[j-1][i])
	 		+F2*(u[j+1][i]-u[j-2][i]);
	 		txz[j][i] = txz[j][i] +
	 		c55[j][i]*(dtz*duz + dtx*dwx);
	 	  }
	 	}
}

void update_stress_ani(int ifx, int ilx, int jfz, int jlz,
	float **u, float **w, float **txx, float **tzz, float **txz,
	float **c11, float **c55, float **c33, float **c13, 
	float **c15, float **c35, float dtx, float dtz)
{
	int i,j;
	float dux,duz,dwz,dwx;

 	for (j=jfz; j<=jlz; ++j) {
		for (i=ifx; i<=ilx; ++i) {

			dwx = F1*(w[j][i+1]-w[j][i])
				+F2*(w[j][i+2]-w[j][i-1]);

			duz = F1*(u[j][i]-u[j-1][i])
	 					+F2*(u[j+1][i]-u[j-2][i]);

	 				dwz = F1*(w[j+1][i]-w[j][i])
	 					+F2*(w[j+2][i]-w[j-1][i]);

					dux = F1*(u[j][i]-u[j][i-1])
			 		+F2*(u[j][i+1]-u[j][i-2]);

			txx[j][i] = txx[j][i] 
					+ c11[j][i]*dtx*dux 
					+ c13[j][i]*dtz*dwz 
					+ c15[j][i]*(dtz*duz + dtx*dwx);

					tzz[j][i] = tzz[j][i] 
					+ c33[j][i]*dtz*dwz 
					+ c13[j][i]*dtx*dux 
					+ c35[j][i]*(dtz*duz + dtx*dwx); 

			txz[j][i] = txz[j][i] 
					+ c15[j][i]*dtx*dux 
					+ c35[j][i]*dtz*dwz 
					+ c55[j][i]*(dtz*duz + dtx*dwx);
	 	  	}
		i=ilx+1;

		dwx = F1*(w[j][i+1]-w[j][i])
		 		+ F2*(w[j][i+2]-w[j][i-1]);

		duz = F1*(u[j][i]-u[j-1][i])
	 				+ F2*(u[j+1][i]-u[j-2][i]);

	 			dwz = F1*(w[j+1][i]-w[j][i])
			+ F2*(w[j+2][i]-w[j-1][i]);

	 			dux = F1*(u[j][i]-u[j][i-1])
			+ F2*(u[j][i+1]-u[j][i-2]);

		txx[j][i] = txx[j][i] 
			+ c11[j][i]*dtx*dux 
			+ c13[j][i]*dtz*dwz 
			+ c15[j][i]*(dtz*duz + dtx*dwx);

	 			tzz[j][i] = tzz[j][i] 
			+ c33[j][i]*dtz*dwz 
			+ c13[j][i]*dtx*dux 
			+ c35[j][i]*(dtz*duz + dtx*dwx); 
	}
}


void add_p_source(float **txx, float **tzz, float amp, int i, int j)
{
	txx[j][i] = txx[j][i] + amp;
	tzz[j][i] = tzz[j][i] + amp;
}

void add_v_source(float **u, float **w, float amp, int i, int j)
{
	w[j][i] = w[j][i] + amp*0.000001;
}


void add_pw_source_V(float **u, float **w, float **txx, float **tzz, 
	float **txz, float **c11, float **c55, float *a, float sang,
	int il, int ir, int jt, int jb, int ns, int k, float dx,
	float dz, float dt, float *vl, float *vb, float *vr,
	float *ttl, float *ttb, float *ttr)
{
	int i,j;
	float sa, ca, ak, p;
	float tw=0, tu=0, tiw, tiu, tt;
	int ksw, ksu;
	double pi=D_PI;

	/* sang = direction of wave propagation relative to vertical */
	tt=k*dt;

	/* introduce plane wave on left boundary */
	p=sin(sang*pi/180.)/vl[jb];
	for (j=jb+2; j>=jt; --j) {
		sa=sin(asin(p*vl[j]));
		ca=cos(asin(p*vl[j]));

		for (i=il-2; i<il; ++i) {
			if(sang>=0&&sang<90) { 
				tw=tt-ttl[j]-(dx*(i-il+0.75)*sa 
					+ dz*(-0.25)*ca)/vl[j];
				tu=tt-ttl[j]-(dx*(i-il+1.25)*sa 
					+ dz*(-0.75)*ca)/vl[j];
	   		}
	   		if(sang>-90&&sang<0) { 
					tw=tt-ttl[j]-(dx*(i-ir-1.75)*sa 
					+ dz*(-0.25)*ca)/vl[j];
					tu=tt-ttl[j]-(dx*(i-ir-1.25)*sa 
					+ dz*(-0.75)*ca)/vl[j];
	   		}
	   		ksw=floor(tw/dt);
			tiw=fmod(tw,dt)/dt;
	   		ksu=floor(tu/dt);
			tiu=fmod(tu,dt)/dt;
				if (ksw>=0 && ksw <ns-1) {
					ak=(1-tiw)*a[ksw]+tiw*a[ksw+1];
				w[j][i]=-vl[j]*ak*ca;
			} 

			if (ksu>=0 && ksu <ns-1) {
				ak=(1-tiu)*a[ksu]+tiu*a[ksu+1];
				u[j][i]=vl[j]*ak*sa;
			}
		}
	}

	/* introduce plane wave on right boundary */
	p=sin(sang*pi/180.)/vr[jb];
	for (j=jb+2; j>=jt; --j) {
		sa=sin(asin(p*vr[j]));
		ca=cos(asin(p*vr[j]));

		for (i=ir+1; i<ir+3; ++i) {
			if(sang>=0&&sang<90) { 
				tw=tt-ttr[j]-(dx*(i-il+1.75)*sa 
					+ dz*(-0.25)*ca)/vr[j];
				tu=tt-ttr[j]-(dx*(i-il+1.25)*sa 
					+ dz*(-0.75)*ca)/vr[j];
			}
			if(sang>-90&&sang<0) { 
				tw=tt-ttr[j]-(dx*(i-ir-0.75)*sa 
					+ dz*(-0.25)*ca)/vr[j];
				tu=tt-ttr[j]-(dx*(i-ir-1.25)*sa 
					+ dz*(-0.75)*ca)/vr[j];
		   	}
			ksw=floor(tw/dt);
			tiw=fmod(tw,dt)/dt;
			ksu=floor(tu/dt);
			tiu=fmod(tu,dt)/dt;

			if (ksw>=0 && ksw <ns-1) {
				ak=(1-tiw)*a[ksw]+tiw*a[ksw+1];
				w[j][i+1]=-vr[j]*ak*ca;
			}

			if (ksu>=0 && ksu <ns-1) {
				ak=(1-tiu)*a[ksu]+tiu*a[ksu+1];
				u[j][i]=vr[j]*ak*sa;
			}
		}
	}
	 
	/* introduce plane wave on bottom boundary */
	 if(sang>=0&&sang<90) p=sin(sang*pi/180.)/vb[il];
	 if(sang>-90&&sang<0) p=sin(sang*pi/180.)/vb[ir];
	 for (i=il; i<ir+2; ++i){
		sa=sin(asin(p*vb[i]));
		ca=cos(asin(p*vb[i]));

	  	for (j=jb-1; j<jb+3; ++j) {
			if(sang>=0&&sang<90) { 
					tw=tt-ttb[i]-(dx*(+0.75)*sa 
					+ dz*(jb-j-0.25)*ca)/vb[i];
					tu=tt-ttb[i]-(dx*(+1.25)*sa 
					+ dz*(jb-j-0.75)*ca)/vb[i];
	   		}
	   		if(sang>-90&&sang<0) { 
					tw=tt-ttb[i]-(dx*(-1.75)*sa 
					+ dz*(jb-j-0.25)*ca)/vb[i];
					tu=tt-ttb[i]-(dx*(-1.25)*sa 
					+ dz*(jb-j-0.75)*ca)/vb[i];
	   		}
	   		ksw=floor(tw/dt);
			tiw=fmod(tw,dt)/dt;
			ksu=floor(tu/dt);
			tiu=fmod(tu,dt)/dt;

			if (ksw>=0 && ksw <ns-1) {
					ak=(1-tiw)*a[ksw]+tiw*a[ksw+1];
				w[j][i]=-vb[i]*ak*ca;
	   		}
			if (ksu>=0 && ksu <ns-1) {
				ak=(1-tiu)*a[ksu]+tiu*a[ksu+1];
				u[j][i]=vb[i]*ak*sa;
	   		}
		}	
	}

}

void add_pw_source_S(float **u, float **w, float **txx, float **tzz,
	float **txz, float **c11, float **c55, float *a, float sang,
	int il, int ir, int jt, int jb, int ns, int k, float dx,
	float dz, float dt, float *vl, float *vb, float *vr,
	float *ttl, float *ttb, float *ttr)
{
	int i,j;
	float sa, ca, ak, p;
	float tw=0, tu=0, tiw, tiu, tt;
	int ksw, ksu;
	double pi=D_PI;

	/* sang = direction of wave propagation relative to vertical */
	tt=k*dt;

	/* introduce plane wave on left boundary */
	p=sin(sang*pi/180.)/vl[jb];
	for (j=jb+2; j>=jt; --j) {
		sa=sin(asin(p*vl[j]));
		ca=cos(asin(p*vl[j]));

		for (i=il-2; i<il; ++i) {
			if(sang>=0&&sang<90) { 
				tw=tt-ttl[j]-(dx*(i-il+0.75)*sa 
					+ dz*(-0.75)*ca)/vl[j];
				tu=tt-ttl[j]-(dx*(i-il+1.25)*sa 
					+ dz*(-0.25)*ca)/vl[j];
	   		}
			if(sang>-90&&sang<0) { 
				tw=tt-ttl[j]-(dx*(i-ir-1.75)*sa 
					+ dz*(-0.75)*ca)/vl[j];
				tu=tt-ttl[j]-(dx*(i-ir-1.25)*sa 
					+ dz*(-0.25)*ca)/vl[j];
	   		}
	   		ksw=floor(tw/dt);
			tiw=fmod(tw,dt)/dt;
			ksu=floor(tu/dt);
			tiu=fmod(tu,dt)/dt;

			if (ksw>=0 && ksw <ns-1) {
				ak=(1-tiw)*a[ksw]+tiw*a[ksw+1];
				txx[j][i]=(c11[j][i]*sa*sa 
					   +(c11[j][i]-2*c55[j][i])*ca*ca)*ak;
				tzz[j][i]=(c11[j][i]*ca*ca
					   +(c11[j][i]-2*c55[j][i])*sa*sa)*ak;
	   		}
			if (ksu>=0 && ksu <ns-1) {
				ak=(1-tiu)*a[ksu]+tiu*a[ksu+1];
				txz[j][i]=c55[j][i]*(2*ca*sa)*ak;
	   		}
		}	
	}

	/* introduce plane wave on right boundary */
	 p=sin(sang*pi/180.)/vr[jb];

	for (j=jb+2; j>=jt; --j) {
		sa=sin(asin(p*vr[j]));
		ca=cos(asin(p*vr[j]));
		for (i=ir+1; i<ir+3; ++i) {
			if(sang>=0&&sang<90) { 
				tw=tt-ttr[j]-(dx*(i-il+1.75)*sa 
					+ dz*(-0.75)*ca)/vr[j];
				tu=tt-ttr[j]-(dx*(i-il+1.25)*sa 
					+ dz*(-0.25)*ca)/vr[j];
	   		}
			if(sang>-90&&sang<0) { 
				tw=tt-ttr[j]-(dx*(i-ir-0.75)*sa 
					+ dz*(-0.75)*ca)/vr[j];
				tu=tt-ttr[j]-(dx*(i-ir-1.25)*sa 
					+ dz*(-0.25)*ca)/vr[j];
	   		}
			ksw=floor(tw/dt);
			tiw=fmod(tw,dt)/dt;

			ksu=floor(tu/dt);
			tiu=fmod(tu,dt)/dt;

			if (ksw>=0 && ksw <ns-1) {
				ak=(1-tiw)*a[ksw]+tiw*a[ksw+1];
					txx[j][i+1]=(c11[j][i]*sa*sa
							+(c11[j][i]
							-2*c55[j][i])*ca*ca)*ak;
					tzz[j][i+1]=(c11[j][i]*ca*ca
							+(c11[j][i]
							-2*c55[j][i])*sa*sa)*ak;
			} else  {
				txx[j][i+1]=0;
				tzz[j][i+1]=0;
			}

			if (ksu>=0 && ksu <ns-1) {
				ak=(1-tiu)*a[ksu]+tiu*a[ksu+1];
				txz[j][i]=c55[j][i]*(2*ca*sa)*ak;
			} else {
				txz[j][i]=0;
			}
		}
	}
	 
	/* introduce plane wave on bottom boundary */
	p=sin(sang*pi/180.)/vb[il];
	for (i=il; i<ir+2; ++i){
		sa=sin(asin(p*vb[i])); ca=cos(asin(p*vb[i]));
		for (j=jb-1; j<jb+3; ++j) {
			if(sang>=0&&sang<90) { 
				tw=tt-ttb[i]-(dx*(+0.75)*sa 
					+ dz*(jb-j-0.75)*ca)/vb[i];
				tu=tt-ttb[i]-(dx*(+1.25)*sa 
					+ dz*(jb-j-0.25)*ca)/vb[i];
			}
			if(sang>-90&&sang<0) { 
				tw=tt-ttb[i]-(dx*(-1.75)*sa 
					+ dz*(jb-j-0.75)*ca)/vb[i];
				tu=tt-ttb[i]-(dx*(-1.25)*sa 
					+ dz*(jb-j-0.25)*ca)/vb[i];
	   		}
			ksw=floor(tw/dt);
			tiw=fmod(tw,dt)/dt;
			ksu=floor(tu/dt);
			tiu=fmod(tu,dt)/dt;
			if (ksw>=0 && ksw <ns-1) {
				ak=(1-tiw)*a[ksw]+tiw*a[ksw+1];
				txx[j][i]=(c11[j][i]*sa*sa
						+(c11[j][i]
							-2*c55[j][i])*ca*ca)*ak;
				tzz[j][i]=(c11[j][i]*ca*ca
						+(c11[j][i]
							-2*c55[j][i])*sa*sa)*ak;
	   		}

			if (ksu>=0 && ksu <ns-1) {
				ak=(1-tiu)*a[ksu]+tiu*a[ksu+1];
				txz[j][i]=c55[j][i]*(2*ca*sa)*ak;
	   		}
		}	
	}

}

void make_snap(float **u, float **w, float sx, float sz, 
		float dx, float dz, int nxpadded, int nzpadded,
		int nx, int nz, int k, float t, float fx, float fz,
		segy tr, FILE *sneisfp, int *wbc)
{
	
	/* SEGY fields */
	long tracl=0;		/* trace number within a line */
	long tracr;		/* trace number within a reel */
	int i,j;
	
	tr.sx = sx;
	tr.selev = -sz;
	tr.trid = 1;
	tr.fldr = k;
	tr.ns = nz;
	tr.dt = 1000 * dz;
	tr.d2 = dx ;
	tr.f1 = fz ;
	tr.f2 = fx ;
	tr.delrt = t; /* store time of snapshot */

	if(&tracl==NULL) tracl = 0;
	tracr = 0;
	tr.trid = TINLIN; /* inline particle velocity component*/
	for (i=0 ; i < nxpadded-wbc[3]-wbc[1] ; ++i){
		++tracl;
		++tracr;			

		tr.gx = fx+dx*(i-wbc[1]);
		tr.tracl = tracl;
		tr.tracr = tracr;
		tr.offset = tr.gx-tr.sx ;

		for(j=0; j<nzpadded-wbc[2]-wbc[0]; ++j){
			tr.data[j] = u[j+wbc[0]][i+wbc[1]];
		}
		fputtr(sneisfp , &tr);
	}
	tr.trid = TVERT ; /* vertical particle velocity component */
	for (i=0 ; i < nxpadded-wbc[3]-wbc[1] ; ++i){
		++tracl;
		++tracr;

		tr.gx = fx+dx*(i-wbc[1]);
		tr.offset = tr.gx-tr.sx ;
		tr.tracl = tracl;
		tr.tracr = tracr;

		for(j=0; j<nzpadded-wbc[2]-wbc[0]; ++j){
			tr.data[j] = w[j+wbc[0]][i+wbc[1]];
		}
		fputtr(sneisfp , &tr);
	}
}

void make_stress_snap(float **txx, float **tzz, float **txz, 
	float sx, float sz, 
	float dx, float dz, int nxpadded, int nzpadded, int k, float t,
	float fx, float fz,
	segy tr, FILE *sneisfp, int *wbc)
{
	
	/* SEGY fields */
	long tracl=0;		/* trace number within a line */
	long tracr=0;		/* trace number within a reel */
	int i,j;
	
	tr.sx = sx;
	tr.selev = -sz;
	tr.trid = 1;
	tr.fldr = k;
	tr.ns = nzpadded ;
	tr.dt = 1000 * dz;
	tr.d2 = dx ;
	tr.f1 = fz ;
	tr.f2 = fx ;
	tr.delrt = t ; /* store time of snapshot */

	if(&tracl==NULL) tracl = 0;
	tracr = 0;
	tr.trid = TINLIN ; /*txx-component in line component*/
	for (i=0 ; i < nzpadded-wbc[3] ; ++i){
		++tracl;
		++tracr;			

		tr.gx = fx+dx*(i-wbc[1]);
		tr.tracl = tracl;
		tr.tracr = tracr;
		tr.offset = tr.gx-tr.sx ;

		for(j=0; j<nxpadded-wbc[2]; ++j){
			tr.data[j] = txx[j+wbc[0]][i+wbc[1]];
		}
		fputtr(sneisfp , &tr);
	}
	tr.trid = TVERT ; /*tzz-component*/
	for (i=0 ; i < nzpadded-wbc[3] ; ++i){
		++tracl;
		++tracr;			

		tr.gx = fx+dx*(i-wbc[1]);
		tr.offset = tr.gx-tr.sx ;
		tr.tracl = tracl; tr.tracr = tracr;

		for(j=0; j<nxpadded-wbc[2]; ++j){
			tr.data[j] = tzz[j+wbc[0]][i+wbc[1]];
		}
		fputtr(sneisfp , &tr);
	}
	tr.duse = TXLIN ; /* txz-component using crossline label */
	for (i=0 ; i < nzpadded-wbc[3] ; ++i){
		++tracl;
		++tracr;			

		tr.gx = fx+dx*(i-wbc[1]);
		tr.tracl = tracl;
		tr.tracr = tracr;
		tr.offset = tr.gx-tr.sx ;

		for(j=0; j<nxpadded-wbc[2]; ++j){
			tr.data[j] = txz[j+wbc[0]][i+wbc[1]];
		}
		fputtr(sneisfp , &tr);
	}
}


void make_seis(float **ut, float **wt, float sx, float sz, 
	float dn, float dt, int nn, int nt, float dd,
	 float fx, float fz,
	segy tr, FILE *fp, int b1, int b2, int ivh, int verbose, float favg)
{
	
	/* SEGY fields */
	long tracl=0;		/* trace number within a line */
	long tracr=0;		/* trace number within a reel */
	int i,j,k;
	

	tr.sdepth = 0;
	tr.sx = sx;
	tr.selev = -sz ;
	tr.ns = nt ;
	tr.d2 = dn ;
	tr.dt = dt * 1000000.0 ;
	tr.ntr = 2*(nn-b2-b1);
	tr.delrt = -NINT(1000.0/(2*favg)) ;
	
	if (verbose) warn("%d %d %d %d %f",tr.ns,tr.tracl,tr.tracr,tr.ntr);

	if(ivh==1) {  /* surface seismic */
		tr.gelev = -dd ;
		tr.f2 = fx ;
	}
	if(ivh==2) {  /* VSP seismic */
	  tr.gx = dd ;
	  tr.f2 = fz ;
	}
	
	tracl = tracr = 0; k=0;
	tr.duse = 2 ; /*u-component*/
	for (i=b1 ; i < nn-b2 ; ++i){
		++tracl; ++tracr;			
		tr.tracl = tracl; tr.tracr = tracr; tr.tracf=i;
		if(ivh==1) {  /* surface seismic */
		  tr.gx = fx+dn*k ;
		  tr.offset = tr.gx-tr.sx ;
		}
		if(ivh==2) {  /* VSP seismic */
		  tr.gelev = fz-dn*k ;
		  tr.offset = tr.sx-tr.gelev ;
		}
		 for(j=0; j<nt; ++j){tr.data[j] = ut[j][i];}
		fputtr(fp , &tr); k++;
	}
	k=0;
	tr.duse = 1 ; /*w-component*/
	for (i=b1 ; i < nn-b2 ; ++i){
		++tracl; ++tracr;			
		tr.tracl = tracl; tr.tracr = tracr; tr.tracf=i;

		if(ivh==1) {  /* surface seismic */
			tr.gx = fx+dn*k ;
			tr.offset = tr.gx-tr.sx ;
		}
		if(ivh==2) {  /* VSP seismic */
			tr.gelev = fz-dn*k ;
			tr.offset = tr.sx-tr.gelev ;
		}
		for(j=0; j<nt; ++j){
			tr.data[j] = wt[j][i];
		}

		fputtr(fp , &tr);
		k++;
	}
}


/* absorbing bc for top*/	
void abs_bc_top(float **u, float **w, float **txx, float **tzz,
			float **txz, int nxpadded, int nzpadded, float *r,
			int *bc, int *wbc, int ifx, int ilx)
{
	int i,j; float fac;
		
	for (j=0; j<wbc[0]; ++j){
		fac=r[(wbc[0]-j-1)];
		for (i=ifx; i<ilx; ++i){
				u[j][i]*=fac;
			w[j][i]*=fac;
				txx[j][i]*=fac;
			tzz[j][i]*=fac;
			txz[j][i]*=fac;
		}
	}
}

/* absorbing bc for left side*/	
void abs_bc_left(float **u, float **w, float **txx, float **tzz,
			float **txz, int nxpadded, int nzpadded, float *r,
			int *bc, int *wbc, int jfz, int jlz)
{
	int i,j;
	float fac;
		
	for (j=jfz; j<jlz; ++j){
		for (i=0; i<wbc[1]; ++i){
			fac=r[(wbc[1]-i-1)];
				u[j][i]*=fac;
			w[j][i]*=fac;
				txx[j][i]*=fac;
			tzz[j][i]*=fac;
			txz[j][i]*=fac;
		}
	}
}

/* absorbing bc for bottom*/	
void abs_bc_bot(float **u, float **w, float **txx, float **tzz,
			float **txz, int nxpadded, int nzpadded, float *r,
			int *bc, int *wbc, int ifx, int ilx)
{
	int i,j;
	float fac;
		
	for (j=nzpadded-wbc[2]; j<nzpadded; ++j){
		fac=r[(j-nzpadded+wbc[2])];
	  	for (i=ifx; i<ilx; ++i){
				u[j][i]*=fac;
			w[j][i]*=fac;
				txx[j][i]*=fac;
			tzz[j][i]*=fac;
			txz[j][i]*=fac;
		}
	}		
}

/* absorbing bc for right side*/	
void abs_bc_right(float **u, float **w, float **txx, float **tzz,
			float **txz, int nxpadded, int nzpadded, float *r, 
			int *bc, int *wbc, int jfz, int jlz)
{
	int i,j;
	float fac;
		
	for (j=jfz; j<jlz; ++j){
		for (i=nxpadded-wbc[3]; i<nxpadded; ++i){
				fac=r[(i-nxpadded+wbc[3])];
				u[j][i]*=fac;
			w[j][i+1]*=fac;
				txx[j][i+1]*=fac;
			tzz[j][i+1]*=fac;
			txz[j][i]*=fac;
		}
	}
}
	
/* Free surface boundary condition where stresses are set to zero
   for the order 4 FD operater on a staggered grid */
void fs4s_bc_top(float **txx, float **tzz, float **txz, int nxpadded, int nzpadded)
{
	int i;
		
	for (i=2; i<nxpadded-2; ++i){
	 		txz[1][i]=0.; txz[0][i]=0.; txz[2][i]=0.;
	}
	for (i=2; i<nxpadded-1; ++i){
	 		txx[1][i]=0.; txx[0][i]=0.; 
	 		tzz[1][i]=-tzz[2][i]; tzz[0][i]=-tzz[3][i];
	}
}

/* Free surface boundary condition where velocities satisfy stress free
  surface for the order 4 FD operater on a staggered grid

** requires anisotropy symmetry parallel to free surface in 
** current implementation */

void fs4v_bc_top(float **u, float **w, float **c11,
			float **c55, int nxpadded, int nzpadded)
{
	int i;
	
	/*  	010514-split loop into two four parts as it was 
		in original FORTRAN CODE, previous versions 
		were calculating all vaulues at once	 */

	for (i=2; i<nxpadded-2; ++i){
		u[0][i]=u[2][i]
			+ 0.5*( F1*(w[2][i+1]-w[2][i]) 
			+ F2*(w[2][i+2]-w[2][i-1]) );
	}

	for (i=2; i<nxpadded-1; ++i){
		w[0][i]=w[2][i]
			+ 0.5*(1-2*(c55[2][i]/c11[2][i]))
			*( F1*(u[0][i]-u[0][i-1]) 
			+ F2*(u[0][i+1]-u[0][i-2]) );
	}
	for (i=2; i<nxpadded-2; ++i){
		u[1][i]=u[0][i]
			+ 0.5*( F1*(w[2][i+1]-w[2][i])
				+ F2*(w[2][i+2]-w[2][i-1]) );
	}
	for (i=2; i<nxpadded-1; ++i){
		w[1][i]=w[0][i]
			+ 0.5*(1-2*(c55[2][i]/c11[2][i]))
			*( F1*(u[1][i]-u[1][i-1]) 
				+ F2*(u[1][i+1]-u[1][i-2]) );
	}
}

void sym4s_bc_top(float **txx, float **tzz, float **txz, int nxpadded, int nzpadded)
/* Stress symmetry boundary condition for the TOP boundary
   for the order 4 FD operater on a staggered grid */
{
	int i;
		
	for (i=2; i<nxpadded-2; ++i){
	 		txx[1][i]=txx[2][i];  txx[0][i]=txx[3][i]; 
	 		tzz[1][i]=tzz[2][i];  tzz[0][i]=tzz[3][i];
	 		txz[1][i]=-txz[3][i]; txz[0][i]=-txz[4][i];
	}
}

void sym4v_bc_top(float **u, float **w, int nxpadded, int nzpadded)
/* Velocity symmetry boundary condition for the TOP boundary
   for the order 4 FD operater on a staggered grid */
{
	int i;
		
	for (i=2; i<nxpadded-2; ++i){
	 		u[1][i]=u[2][i];  u[0][i]=u[3][i]; 
	 		w[1][i]=-w[3][i]; w[0][i]=-w[4][i];
	}
}

void sym4s_bc_right(float **txx, float **tzz, float **txz, int nxpadded, int nzpadded)
/* Stress symmetry boundary condition for the RHS boundary
   for the order 4 FD operater on a staggered grid */
{
	int j;
		
	for (j=0; j<nzpadded; ++j){
	 		txx[j][nxpadded]=txx[j][nxpadded-4]; txx[j][nxpadded-1]=txx[j][nxpadded-3]; 
	 		tzz[j][nxpadded]=tzz[j][nxpadded-4]; tzz[j][nxpadded-1]=tzz[j][nxpadded-3];
	 		txz[j][nxpadded-1]=-txz[j][nxpadded-4]; txz[j][nxpadded-2]=-txz[j][nxpadded-3];
	}
}

void sym4v_bc_right(float **u, float **w, int nxpadded, int nzpadded)
/* Velocity symmetry boundary condition for the RHS boundary
   for the order 4 FD operater on a staggered grid */
{
	int j;
		
	for (j=0; j<nzpadded; ++j){
	 		u[j][nxpadded-1]=-u[j][nxpadded-4]; u[j][nxpadded-2]=-u[j][nxpadded-3]; 
	 		w[j][nxpadded]=w[j][nxpadded-4]; w[j][nxpadded-1]=w[j][nxpadded-3];
	}
}

void sym4s_bc_left(float **txx, float **tzz, float **txz, int nxpadded, int nzpadded)
/* Stress symmetry boundary condition for the LHS boundary
   for the order 4 FD operater on a staggered grid */
{
	int j;
		
	for (j=0; j<nzpadded; ++j){
	 		txx[j][1]=txx[j][3]; txx[j][0]=txx[j][4]; 
	 		tzz[j][1]=tzz[j][3];  tzz[j][0]=tzz[j][4];
	 		txz[j][1]=-txz[j][2]; txz[j][0]=-txz[j][3];
	}
}

void sym4v_bc_left(float **u, float **w, int nxpadded, int nzpadded)
/* Velocity symmetry boundary condition for the LHS boundary
   for the order 4 FD operater on a staggered grid */
{
	int j;
		
	for (j=0; j<nzpadded; ++j){
	 		u[j][1]=-u[j][2]; u[j][0]=-u[j][3]; 
	 		w[j][1]=w[j][3]; w[j][0]=w[j][4];
	}
}

void sym4s_bc_bot(float **txx, float **tzz, float **txz, int nxpadded, int nzpadded)
/* Stress symmetry boundary condition for the BOTTOM boundary
   for the order 4 FD operater on a staggered grid */
{
	int i;
		
	for (i=2; i<nxpadded-2; ++i){
/*	 		txx[nzpadded-2][i]=txx[nzpadded-4][i];  txx[nzpadded-1][i]=txx[nzpadded-5][i];  */
/*	 		tzz[nzpadded-2][i]=-tzz[nzpadded-4][i]; tzz[nzpadded-1][i]=-tzz[nzpadded-5][i]; */
/*	 		txz[nzpadded-2][i]=-txz[nzpadded-3][i]; txz[nzpadded-1][i]=-txz[nzpadded-4][i];} */
	 		txx[nzpadded-2][i]=txx[nzpadded-4][i];  txx[nzpadded-1][i]=txx[nzpadded-5][i]; 
	 		tzz[nzpadded-2][i]=tzz[nzpadded-4][i]; tzz[nzpadded-1][i]=tzz[nzpadded-5][i];
	 		txz[nzpadded-2][i]=-txz[nzpadded-3][i]; txz[nzpadded-1][i]=-txz[nzpadded-4][i];
	}
}

void sym4v_bc_bot(float **u, float **w, int nxpadded, int nzpadded)
/* Velocity symmetry boundary condition for the BOTTOM boundary
   for the order 4 FD operater on a staggered grid */
{
	int i;
		
	for (i=2; i<nxpadded-2; ++i){
/*	 		u[nzpadded-2][i]=u[nzpadded-4][i]; u[nzpadded-1][i]=u[nzpadded-5][i];  */
/*	 		w[nzpadded-2][i]=-w[nzpadded-3][i];w[nzpadded-1][i]=-w[nzpadded-4][i];} */
	 		u[nzpadded-2][i]=u[nzpadded-4][i]; u[nzpadded-1][i]=u[nzpadded-5][i]; 
	 		w[nzpadded-2][i]=-w[nzpadded-3][i]; w[nzpadded-1][i]=-w[nzpadded-4][i];
	}
}

void calc_area(float vmax, float dtx, float dtz, 
		int *limits, int i, int j, int k, int nxpadded, int nzpadded)
{
	 limits[0] = i - NINT(vmax*k*dtx) - 10;
	 limits[1] = i + NINT(vmax*k*dtx) + 10;
	 limits[2] = j - NINT(vmax*k*dtz) - 10;
	 limits[3] = j + NINT(vmax*k*dtz) + 10;

	 if (limits[0] < 2) limits[0] = 2;

	 if (limits[1] > nxpadded-3) limits[1] = nxpadded-3;

	 if (limits[2] < 2) limits[2] = 2;

	 if (limits[3] > nzpadded-3) limits[3] = nzpadded-3;	 
}



void write_chd(int nxpadded, int nzpadded, int nt, float dx, float dz, float dt,
	float fx, float xmax, float fz, float zmax,
	float sx, float sz, float favg, float ts, char *wtype, char *stype,
	int *bc, int qsw, int aniso, float xz,
	FILE *fp, int comp, char *mfile)
{
	int k=1;
	char s[80];
	k=fwrite_chd(fp,"SUEA2DF",k); 
	k=fwrite_chd(fp,"Department of Earth Sciences",k);
	k=fwrite_chd(fp,"Uppsala University",k);
	k=fwrite_chd(fp,"Date:",k);
	k=fwrite_chd(fp," ",k);

	if(comp==1) {
		sprintf(s,"Horizontal line recording at z=%f",xz);
		k=fwrite_chd(fp,s,k);
	}

	if(comp==2) {
		sprintf(s,"Vertical line recording at x=%f",xz);
	  	k=fwrite_chd(fp,s,k);
	}
	sprintf(s,"nxpadded=%d nzpadded=%d nt=%d dx=%f dz=%f dt=%f",nxpadded,nzpadded,nt,dx,dz,dt);
	  k=fwrite_chd(fp,s,k);
	sprintf(s,"fx=%f xmax=%f fz=%f zmax=%f",fx,xmax,fz,zmax);
	  k=fwrite_chd(fp,s,k);
	sprintf(s,"sx=%f sz=%f favg=%f ts=%f",sx,sz,favg,ts);
	  k=fwrite_chd(fp,s,k);
	sprintf(s,"wtype=%s stype=%s",wtype,stype);
	  k=fwrite_chd(fp,s,k);
	sprintf(s,"bc=%d,%d,%d,%d",bc[0],bc[1],bc[2],bc[3]);
	  k=fwrite_chd(fp,s,k);
	k=fwrite_chd(fp," ",k);
	k=fwrite_chd(fp,"Model file:",k);

	sprintf(s," %s",mfile);

	k=fwrite_chd(fp,s,k);

	k=fwrite_chd(fp," ",k);

	if(qsw==1) k=fwrite_chd(fp,"Attenuation included",k);

	if(aniso==1) k=fwrite_chd(fp,"Anisotropy included",k);

	while(k<=40) k=fwrite_chd(fp," ",k);
}	
	
int fwrite_chd(FILE *fp, char *s, int k)
{
	int i,n;	

	n=fprintf(fp,"C %s",s);

	for(i=1; i<80-n-1; ++i) fputs(" ",fp); 

	fprintf(fp,"%s\n"," ");

	return(k+1);
}


void write_grid(float **u, float **w, float **txx, float **tzz,
			float **txz, int is, int js)
{
	int i,j;
		
	for (i=is-2; i<is+3; ++i){
		fprintf(stderr,"%d ",i);
	}

	fprintf(stderr,"\n");

	for (j=js-2; j<js+3; ++j){
	   	for (i=is-2; i<is+3; ++i){
			fprintf(stderr,"%e ",u[j][i]);
		}
		fprintf(stderr,"%d \n",j);
	}

	fprintf(stderr,"w \n");

	for (j=js-2; j<js+3; ++j){
		for (i=is-2; i<is+3; ++i){
			fprintf(stderr,"%e ",w[j][i]);
		}
		fprintf(stderr,"%d \n",j);
	}

	fprintf(stderr,"txx \n");

	for (j=js-2; j<js+3; ++j){
		for (i=is-2; i<is+3; ++i){
			fprintf(stderr,"%e ",txx[j][i]);
		}
		fprintf(stderr,"%d \n",j);
	}

	fprintf(stderr,"tzz \n");

	for (j=js-2; j<js+3; ++j){
		for (i=is-2; i<is+3; ++i){
			fprintf(stderr,"%e ",tzz[j][i]);
		}
		fprintf(stderr,"%d \n",j);
	}

	fprintf(stderr,"txz \n");

	for (j=js-2; j<js+3; ++j){
		for (i=is-2; i<is+3; ++i){
			fprintf(stderr,"%e ",txz[j][i]);
		}
		fprintf(stderr,"%d \n",j);
	}

}

void pad_float_2D(int n1, int n2, int *padval, float **in, float **out) 
{
/*************************************************************************
pad_float_2D - padd a 2D array of floats
**************************************************************************

Input:
n1		size of the 1 (fast) dimension
n2		size of the 2 (slow) dimension
padval[4]	paddings  	padval[0] beginning pad in dimension 1 
				padval[1] beginning pad in dimension 2
				padval[2] ending pad in dimension 1
				padval[3] ending pad in dimension 2
in[][]		n1 by n2 array of input floats
Output:
out[][]		n1padded by n2padded array of input floats
		where the padded sizes are given by
		n1padded = n1 + padval[1] + padval[3]
		n2padded = n2 + padval[0] + padval[2]
**************************************************************************
Notes:
**************************************************************************
Author: CWP: John Stockwell, 2005
**************************************************************************/

	int ix1;	/* counter in dimension 1 */
	int ix2;	/* counter in dimension 2 */
	int n1padded;	/* size of dimension 1 for padded array */
	int n2padded;	/* size of dimension 2 for padded array */

	/* calculate padded six2es */
	n1padded = n1 + padval[1] + padval[3];
	n2padded = n2 + padval[0] + padval[2];

	/* zero out the out[][] array */
	memset((void *) out[0], 0, n1padded*n2padded*FSIZE);

	/* copy in[][] into out[][]  */
	for (ix2=padval[0]; ix2 < n2; ++ix2) 
	  	for (ix1=padval[1]; ix1<n1; ++ix1) 
			out[ix2][ix1]=in[ix2-padval[0]][ix1-padval[1]];

	/* pad upper left corner */
	for (ix2=0; ix2 < padval[0]; ++ix2) 
	  	for (ix1=0; ix1< padval[1]; ++ix1) 
			out[ix2][ix1]=in[0][0];

	/* pad upper side */
	for (ix2=padval[0]; ix2 < n2; ++ix2) 
	  	for (ix1=0; ix1< padval[1]; ++ix1) 
			out[ix2][ix1]=in[ix2-padval[0]][0];

	/* pad upper right corner */
	for (ix2=0; ix2 < padval[0]+padval[2]; ++ix2) 
	  	for (ix1=0; ix1< padval[1]; ++ix1) 
			out[ix2+n2][ix1]=in[n2-1][0];

	/* pad left side */
	for (ix2=0; ix2 < padval[0]; ++ix2) 
	  	for (ix1=padval[1]; ix1< n1 ; ++ix1) 
			out[ix2][ix1]=in[ix2][ix1-padval[1]];

	/* lower left corner */
	for (ix2=0; ix2 < padval[0]; ++ix2) 
	  	for (ix1=n1; ix1< n1+padval[1]+padval[3] ; ++ix1) 
			out[ix2][ix1]=in[n2-1][n1-1];

	/* pad lower side */
	for (ix2=padval[0]; ix2 < n2; ++ix2) 
	  	for (ix1=n1; ix1< n1+padval[1]+padval[3]; ++ix1) 
			out[ix2][ix1]=in[ix2-padval[0]][n1-1];

	/* pad right side */
	for (ix2=n2; ix2 < n2+padval[0]+padval[2]; ++ix2) 
	  	for (ix1=padval[1]; ix1< n1 ; ++ix1) 
			out[ix2][ix1]=in[n2-1][ix1-padval[1]];

	/* pad lower right corner */
	for (ix2=n2; ix2 < n2+padval[0]+padval[2]; ++ix2) 
	  	for (ix1=n1; ix1< n1+padval[1]+padval[3] ; ++ix1) 
			out[ix2][ix1]=in[n2-1][n1-1];

}
